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Abstract 



The three-dimensional SU(3) spin model is an effective Polyakov loop 
model for QCD at nonzero temperature and density. It suffers from a sign 
problem at nonzero chemical potential. We revisit this model using com- 
plex Langevin dynamics and assess in particular the justification of this 
approach, using analyticity at small /j? and the criteria for correctness de- 
veloped recently. Finite-stepsize effects are discussed in some detail and 
a higher-order algorithm is employed to eliminate leading stepsize correc- 
tions. Our results strongly indicate that complex Langevin dynamics is 
reliable in this theory in both phases, including the critical region. This is 
in sharp contrast to the case of the XY model, where correct results were 
obtained in only part of the phase diagram. 
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1 Introduction 



The phase structure of QCD as temperature and baryon chemical potential are 
varied has not yet been determined from first principles [1] . Due to the presence 
of the sign problem (at nonzero chemical potential the fermion determinant is 
complex), the cornerstone of numerical lattice gauge theory, importance sam- 
pling, is not applicable. In the past decade several approaches have been put 
forward to access the phase diagram at small quark chemical potentials /i and 
at temperatures near the transition temperature between the confined and the 
quark-gluon plasma phase. While in agreement when // < T, none of these meth- 
ods can be extended to larger /i values (see Refs. [2l|3] for recent reviews). Given 
that current and upcoming experiments at the Relativistic Heavy Ion Collider at 
BNL, the Large Hadron Collider at CERN and the Facility for Antiproton and 
Ion Research at GSI aim to map out the phase boundaries in the QCD phase 
diagram by colliding heavy ions at relativistic speeds, there is ample motivation 
to study the sign problem and ways to resolve it, both in QCD and in related 
theories. 

There are several methods which allow the sign problem to be eliminated 
altogether, but these are not universally applicable. In some theories it is possible 
to group degrees of freedom together in such a way that the sign problem is 
manifestly absent. This is the idea behind the meron cluster algorithm |1] and it 
has been applied recently to random matrix theory at finite chemical potential [5] . 
A constrained sampling of field space, yielding a joint probability distribution for 
only a small number of observables, is the notion behind the factorization/density 
of states/histogram approaches P-[9]. Sometimes it is possible to reformulate a 
theory in terms of dual variables in a sign-problem free manner [TOl[Tl]. Recent 
successful applications of this have been to models derived from QCD in combined 
strong-coupling and heavy-quark expansions [T2l - [T8] . 

In this paper we consider complex Langevin dynamics, a numerical algorithm 
not relying on importance sampling but instead on a complexification of the 
degrees of freedom, resulting in new ways to explore field space [T9H2T]. We 
consider the three-dimensional SU(3) spin model at nonzero density, an effective 
Polyakov loop model which follows from the QCD Lagrangian in a combined 
strong-coupling and heavy-quark expansion and one of the first QCD-related 
models addressed with complex Langevin dynamics [221I2S] • Our reason to revisit 
this model is partly due to the recent discussion of Gattringer, who showed how a 
reformulation in terms of fluxes eliminates the sign problem [15]. Moreover, given 
that our understanding of complex Langevin dynamics has steadily improved in 
the past years we consider it worthwhile to reconsider the model and 

apply recently developed tools [2111321131] to assess the applicability of complex 
Langevin dynamics in detail, something that was not undertaken in the classic 
papers [221123]. 

This paper is organized as follows. In the next section we introduce the SU(3) 
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model and summarize some basic results at finite density. The complex Langevin 
equations are given in Sec. |3l Besides the standard lowest-order discretization, 
we also describe a higher-order algorithm to eliminate leading stepsize correc- 
tions [35]. In Sec. H] we discuss our current understanding of the applicability 
of complex Langevin dynamics at finite density and review the various ways in 
which the outcome of a complex Langevin process can be assessed, in particular 
when the exact result is not available. Sec. [5] constitutes the main part of the 
paper. Here we present a variety of numerical results assessing the applicability 
of complex Langevin dynamics in this model, both in the disordered and the 
ordered phase. We also demonstrate that the higher-order algorithm eliminates 
most of the stepsize dependence. Conclusions are drawn in Sec. [61 The higher- 
order algorithm is discussed in some more detail in Appendix [Aj while Appendix 
[B] can be used to scrutinize the stepsize dependence and criteria for correctness. 
A brief account of part of this work has appeared in Ref. |36j . 



2 SU(3) spin model 

We consider the three-dimensional SU(3) spin model at nonzero chemical poten- 
tial, with the action [22] 

S = Sb + Sf, (2.1) 

where 

3 

SB = -f3j2Yl UxTt Ul^, + Tr f/jTr U,+o) , (2.2) 

X u = l 

Sp = -hJ2 (e^Tr + e-'^Tr t/]) . (2.3) 

X 

The model can be thought of as an effective dimensionally reduced version of 
QCD, where TtUx represents the trace of the Polyakov loop; the f^'s are SU(3) 
matrices living on a three-dimensional lattice (we use periodic boundary con- 
ditions). The first term then represents the gluon contribution with effective 
coupling /3, while the second term represents heavy quarks, with coupling h. 
Chemical potential favours quarks over anti-quarks, resulting in a complex action, 
Sp{fi) = Spi—f^*)- The fermion term is a simplified version of the contribution 
derived in the heavy dense limit [27]. The partition function. 




is even in /i due to charge conjugation invariance. Here / denotes the free energy 
density and Q is the three-dimensional volume. 

The phase structure of this theory has been studied in Refs. [221122], using 
both complex Langevin dynamics and mean-field theory. Recently it has also been 
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investigated using a reformulation of the theory which is sign-problem free [T6] . 
For small the theory has a disordered (confined) phase for smaller /3 values, and 
an ordered (deconfined) phase for larger /3 values. The two phases are separated 
by a first-order phase transition. This is the case for vanishing and small chemical 
potential. With increasing chemical potential, the transition weakens and turns 
into a crossover at a critical endpoint. For larger /i, there is a crossover only. 

We will also consider two closely related models which have a real action, 
namely the model with imaginary chemical potential, /i = z/ij, and the phase- 
quenched model, obtained by discarding the imaginary part of the action, such 
that 

= -h cosh /X (Tr + Tr Ul) . (2.5) 

X 

In contrast to QCD, the SU(3) spin model does not have a Silver Blaze prob- 
lem. The Silver Blaze problem [37] refers to the region in the phase diagram where 
the chemical potential is nonzero but bulk thermodynamic observables, such as 
the pressure and the density, are /i-independent. This /i-independence requires 
a precise cancelation which can be highly non-trivial in a numerical approach, 
as can be seen from studies of the eigenvalues of the Dirac operator [371[38]. We 
note here that complex Langevin dynamics has been shown to be able to solve the 
Silver Blaze problem, in the relativistic Bose gas [28l|29] and in one-dimensional 
QCD [33]. To see that the Silver Blaze region is absent in this model, consider 
the density, 

in) = = {he^TiU, - he-^TiUl) . (2.6) 

A nonzero density induces a difference between i^iUx) and (Tr?7j). On the 
other hand, in the Silver Blaze region, (n) = and (Trf/^.) = (YiU].). It is clear 
from the expression above that it is not possible to simultaneously satisfy these 
conditions when yU 7^ 0, hence the Silver Blaze region is absent Q Similarly, the 
density in the phase-quenched theory, 

(r2)pq = h sinh /i (Tr + Tr Ul)^^ , (2.7) 

is nonzero as soon as /i > 0. 

The severeness of the sign problem is conventionally estimated via the expec- 
tation value of the complex phase factor e**^ = /\e~^\ in the phase quenched 
theory, 

(e^^)pq = ^ = e-^(^-^-). (2.8) 

-'^For completeness, we recall that l^rUx) and (Trt/J) are both real in the full theory, while 
at imaginary /i the real parts are equal and the imaginary parts are opposite. In the phase- 
quenched theory, they are real and identical. Note also that the fermion contribution breaks 
the Z(3) symmetry of the bosonic sector, hence (TrJXr) and (Trt/J) are never strictly zero. 
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The full and the phase-quenched theory differ as soon as is nonzero, which can 

be seen by performing a Taylor series expansion around /i = 0. To second order 
in /i, the free energy densities read 

/(/i) = /(O) - (ci + C2h) V + 0{^i'), (2.9) 

/pq(/^) = /(0)-CiV + C^(/), (2.10) 

with 

^i = ^E(Trt/.Uo> (2-11) 

X 

^2 = ^ E i^^ - ^3 Tr (t/, - f/t)>^^^ . (2.12) 

Since C2 is neg ative [Tr {U^ - t/J) is imaginary], / — /pq > 0, as it should be. 
Similarly, (n) < (n)pq. 



3 Discretized complex Langevin dynamics 

In order to solve the complex Langevin evolution numerically, we follow Ref. [22] 
and diagonalize the SU(3) matrices in terms of the angles 0i,2, such that 



Trf/, 
TtUI 



Mix 



+ e 



+ e 



(3.1) 
(3.2) 



We then have to include the reduced Haar measure and consider the partition 
function 

-ScS 



n 



hx a<p2x e 



where 

with 

Sh = 



S, 



cS 



Sb -\' Sp -\- S 



sm 



>2x 



sm 



+ 



'2x 



sm 



>\x 



+ 



nx 



(3.3) 
(3.4) 

. (3.5) 



We note that it is also possible to implement complex Langevin dynamics directly 
for the SU(3) matrices, see e.g. Refs. [231127] . 

Langevin dynamics provides a stochastic update for the angles ^ax {fl = 1, 2), 
according to 

70ax — Kax ~\~ ^axi K^x — 771 ; (3-6) 



d(j)a 
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where denotes the Langevin time, Kax is the drift term, and the noise satisfies 

{riax) = 0, iVaxVa'x') = '^^aa'^xx'- (3.7) 

When the action and hence the drift terms are complex, the angles do not 
remain real under the Langevin evolution. We therefore write (pax = 4>^x + ^4>\xi 
and consider the following complex Langevin equations, using real noise, 

- ^-n - -Re ^^""^ i 

Q^^ax — ^^ax ^ 'lax, J^^ax — , V'-'''^ 



■'ax 



Q^(Pax - ^ax, ^ax " ^^1 



(3.9) 



After complexification, we write instead of f/J in the remainder. 

To solve these equations numerically, Langevin time is discretized a.s = en, 
where e is the Langevin time step. The standard algorithm discretizing Eq. 03.61) 
read^ 

(paxin + 1) = (paxin) + eK[(pax{n)] + y/er]ax{n), (3.10) 

where 

iVaxin)) = 0, {r]ax{n)r]a'x'{n')) = 26aa'Sxx'Snn'- (3.11) 

The contribution to the drift term from the Haar measure requires careful inte- 
gration. For this we use the adaptive stepsize algorithm of Ref. [30] . 

It is well-known that Langevin dynamics has finite stepsize corrections, which 
are linear in e in the lowest-order discretization given above [21]. It is therefore 
necessary to extrapolate to zero stepsize. In our previous work [2714M] . we have 
only considered the lowest-order algorithm. However, motivated by the results 
to be presented below, we implemented a higher-order algorithm to improve the 
stepsize dependence. A standard Runge-Kutta scheme, where the drift terms are 
improved but the noise is kept as above, will not remove the leading stepsize 
correction [39]. Instead, it is necessary to modify the noise terms as well. We use 
the algorithm proposed in Ref. [35] for real Langevin dynamics, which is explicit 
and easy to implement H It takes the following form 

i^ax{n) = (pax{n) + ]^eK[<pax{n)\, 

1 3 _ 

i^ax{n) = <pax{n) + -eK[(pax{n)] + -^/eaax{n), 

(Pax{n + 1) = <pax{n) + [K[^ax{n)] + 2K[iJax{n)]^ + V~eaax{n). (3.12) 
Here otax (n) is a random variable taken according to 

1 

aaxin) = -aax{n) + —Ux{n), (3.13) 
2 D 



^Complexification is obvious and we do not give the discretized equations explicitly. 
•^For other approaches, see e.g. Refs. [401141] . 
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while aaxin) and ^axin) are independent Gaussian random variables with variance 
2 and vanishing mean, i.e., 

{aax{n)aa'x'{ri)) = {^ax{n)^a'x'{n')) = 26aa'Sxx'Snn' , 

{aaxin)^a'x'in)) = (a^xH) = (Uin)) = 0. (3.14) 

In Ref. [35] it was shown analytically, for the case of a real drift term, that 
with this update the remaining correction is O(e^) for a system with one degree 
of freedom and (9(e^/^) for a coupled system. In Appendix |X] we discuss this 
algorithm in some more detail. 



4 Justification and criteria for correctness 

In the case of a real action, it can be shown that stochastic quantization and 
Langevin dynamics is equivalent to standard path integral quantization [21]. As 
is well-known, such a general statement is lacking in the case of a complex action 
[2T] . Indeed, it can occur that under complex Langevin evolution expectation 
values converge to a wrong result [25l[26|[32lH2] - H5] . It is therefore important to 
be able to judge the outcome of a complex Langevin process using assessments 
which are general and can be used in a variety of theories, especially when there 
are no known results to compare with. 

The first assessment employs analyticity at small /x^: observables, which are 
even under charge conjugation, should be analytic as a function of /x^ (in a finite 
volume) [MIITTj. Results at positive /i^ can be compared with those where /i^ < 0, 
i.e. at imaginary potential, obtained using real Langevin dynamics or any other 
standard approach. This test is limited to small chemical potentials. 

A more formal justification of the complexified dynamics can be found in 
Refs. [3ll[3l]. Here we summarize that discussion briefiy in order to arrive at the 
criteria for correctness developed in Ref. [34] . We consider expectation values with 
respect to the real and positive probability distribution P[0^, 0^; ^9], sampled by 
the stochastic process. 



Here d is the Langevin time. With the help of the Langevin equations for 0^ and 
0\ one finds the Fokker-Planck equation for 0^; t?]. 



aP[0^0I;^] _ ^, d 



L^P[0^,0i;^] 



-^K\ (4.2) 



We also consider expectation values with respect to a complex weight p[0; 

/Z)0p[0;^]O[0] 
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where is real and p[0; satisfies 



d_ dS_ 
d(j) d(j) 



(4.4) 



Note that the Fokker-Planck operators L'^ , acting on the real density 0^; 
and Lq, acting on the complex density ^9], should be distinguished. Further- 
more, Eq. (14. 4p has a stationary solution, ~ exp(— S"), whereas for Eq. (14. 2 p 
no generic stationary solution is known. 

Employing that the only permissible observables are holomorphic and making 
use of partial integration, one can show that expectation values with respect to 
the two densities are equal, 

(0)p(,) = (0),(^). (4.5) 

If subsequently one can show that p[(f); reaches the unique stationary solution 
~ exp(— S") in the limit that — )■ oo, the use of complex Langevin dynamics is 
justified [3T] . 

The equivalence in Eq. (14. 5 p relies on the ability to do partial integration 
without receiving contributions from boundary terms, i.e. the distributions should 
be well localized and decay strongly. It was shown in Ref. [31] that this condition 
can be expressed as a set of criteria on holomorphic observables O, which take 
the form 

{LO) = 0. (4.6) 
Here L denotes the Langevin operator 

I^k]' K^-'-l, (4,7) 



(90 J (90 

which depends on holomorphic degrees of freedom 0. Although it differs from L 
and Lq, the action of L on holomorphic observables agrees with that of L. The 
expectation value in Eq. (14. 6 p is taken with respect to the weight P in the limit 
that the Langevin process has equilibrated — t- oo). In principle, the criteria 
(14. 6 p should be satisfied for a large enough set of holomorphic observables [3i 
Adapting this to the model at hand, L reads 



d \ d 



We will consider only local observables and denote these as O[0i3;, 02a;] = Ox- We 
can then write 

lOx = (of + KaxO^j) , (4.9) 



where 



ax 
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In terms of the real and imaginary parts, this yields explicitly 

Re LO. =5^ (Re Of + K^ReO^J - KllmO'j) , (4.11) 

a 

Im LO. = J2 Of + Kfjm + K^Re ) , (4.12) 

a 

and the criteria read 

(ReLO^) = 0, (ImLO^) = 0. (4.13) 

5 Results and justification 

In this section we present a number of results obtained with complex Langevin 
dynamics. As mentioned earlier, our goal is not deliver a detailed study of critical 
properties and the phase structure; rather the aim of this study is to assess the 
reliability of the complex Langevin algorithm, using the criteria discussed in the 
previous section. 




Figure 1: (Tr {Ux + U~^) /2) as a function of /3 at /i = and h = 0.02 on a 10'^ 
lattice, using real Langevin dynamics. 

The results we show here are obtained using a relatively small value for the 
fermion coupling, h = 0.02, so that there is a clear transition between the ordered 
and the disordered phase. We consider /3 values between 0.12 and 0.139; the 
critical /3 value for the fermion coupling we use is around 0.1324 at /i = 0. This is 
demonstrated in Fig.Hl where we show (Tr (f/^. + U~^) /2) = (Tr [4.) as a function 
of /3 at /i = on a 10'^ lattice. 
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Figure 2: Analyticity in yU^: (Tr {Ux + U~^) /2) as a function of yU^ for various /3 
values with /i = 0.02 on a 10^ lattice. Data at imaginary /i (with /i^ < 0) has 
been obtained with real Langevin dynamics, data at real fi (with /i^ > 0) with 
complex Langevin dynamics. 

As a first test, we probe the transition by varying /i instead of /3. As mentioned 
above, observables which are invariant under charge conjugation should, in a finite 
volume, be analytic in /i^ ^61117] . This yields the possibility to compare results 
at positive /i^ with those at negative /i^, corresponding to imaginary potential. 
Since in this case the action is real, real Langevin dynamics can be used, which is 
theoretically well founded. In Fig. [2] we show (Tr {Ux + U~^) /2) as a function of 
yU,^ for eight different /3 values. We observe smooth behaviour as /i^ is increased. 
This is an indication that complex Langevin dynamics works well. We note that 
this is true in both phases as well as in the transition region. This is in contrast 
to the case of the XY model recently studied using complex Langevin dynamics, 
where correct results were obtained in only part of the phase diagram [32] • 

The strength of the transition weakens as /i^ increases and, vice versa, in- 
creases as /i^ decreases [T6l[22l[23]. This can be seen in Fig. [3l where we show 
the Langevin time evolution of (Tr {Ux + U~^) /2) at /i^ = -0.65 and /3 = 0.134 
(left) and /x^ = 0.1 and /3 = 0.132 (right). We observe clear first order behaviour 
at /i^ = —0.65, while at fi^ = 0.1 the transition is much weaker. 

Next we consider the density as a function of fi in the full and the phase- 
quenched theory. In Fig. |4] we show the density for chemical potentials up to 
fi = 3.5, at /3 = 0.125. For this f3 value the model is in the disordered phase at 
smaller fi values and in the ordered phase at larger fi values. The densities in 
the full and phase-quenched theories are similar, but not equal. We recall that 
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Figure 3: Langevin time evolution of (Tr {Ux + U~^) /2) in the transition region, 
at imaginary chemical, /i^ = —0.65 and /3 = 0.134 (left) and real chemical poten- 
tial, /i^ = 0.1 and P = 0.132 (right). The other parameters are as above. 
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Figure 4: Density (n) in the full and the phase- quenched theory as a function 
of /i at /3 = 0.125 and /i = 0.02 on a 10^ lattice. The inset shows a close-up of 
the small /i region. The lines are the predicted linear dependence for small /i, 
evaluated at yU = 0. 

there is no Silver Blaze region in this model. This can be seen in the inset, which 
shows a close-up: the density in the full theory is below the one in the phase- 
quenched theory, but it is nonzero (we have verified that there are no visible 
finite-size effects remaining). The lines indicate the expected linear dependence 
of the densities on /x, using the lowest-order Taylor series expansion, see Eqs. fl2.9[ 
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Figure 5: (Trf/i) and (Tr?7~^) as a function of /i in the full theory. The param- 
eters are as in Fig. HI The inset shows a close-up of the small region. The lines 
are the predicted linear dependence for small /x, evaluated at /i = 0. 

Eini), 

(n) = 2(ci + C2/i)/i/i + C»(/i='), (5.1) 

where the coefficients Ci_2 have been defined in Eqs. fl2.1H 12.120 . In the phase- 
quenched theory the term with C2 is absent. We have computed the coefficients 
and find 

ci = 0.1446(21), C2 = -3.534(72). (5.2) 

Using these coefficients yields the straight lines in the inset of Fig. HJ justifying 
the results of complex Langevin dynamics. 

In Fig. [S]we show (Tr t/) and (Trt/^^) as a function of /i in the full theory, 
using the same parameters as in Fig. |H Recall that (Tr f/) and (Trf/~^) are 
both real and that one expects (Trf/) < (Trf/~^), due to the nonzero density. 
At small /i, the linear dependence on yU can again be expressed in terms of the 
coefficients Ci_2 and we find 

(Trf/) =ci + C2/iAi + 0(Ai'), (5.3) 
(Trt/-i) = Ci-C2V + C(/i'). (5.4) 

This yields the straight lines in the inset of Fig. [5l justifying again the results 
of complex Langevin dynamics. In the phase-quenched theory, (Trf/)pq and 
(Tr ?7~^)pq are equal and slightly below (Tr {U + U~^) /2) in the full theory (not 
shown) . 

The average phase factor, indicating the severeness of the sign problem, is 
shown in Fig.[6]for a typical choice of parameters. The lines indicate the behaviour 
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Figure 6: Average phase factor in the phase- quenched theory (e*'^)pq as a function 
of /i, for various volumes at /3 = 0.125 and h = 0.02. The hues indicates the 
expected behaviour using the leading //^ term at small /i. 

expected at small chemical potential, 

(e^^pq = e"^^^ A/ = / - /p, = -c^h'fi' + Oifi'). (5.5) 

As in preceding studies [27]429| l32 | [33] , we have not observed a correlation between 
the severeness of the sign problem and the efficiency of the complex Langevin 
algorithm. We also note that the average phase factor behaves in a non-monotonic 
manner as a function of fi in the transition region. 

In order to assess complex Langevin dynamics in detail for larger n values, 
we now focus on two points in the phase diagram: /3 = 0.125, /x = 1 in the 
disordered phase and (3 = 0.125, /i = 3 in the ordered phase. To control the 
statistical error we have carried out simulations using a total Langevin trajectory 
of length 10,000 (after discarding the thermalization stage) in the disordered 
phase; in the ordered phase fluctuations are smaller and a Langevin trajectory of 
5,000 is sufficient. Errors are determined with a jackknife analysis. We have used 
a number of stepsizes, from e = 0.001 down to e = 0.00005, employing both the 
standard lowest-order algorithm and the improved higher-order algorithm. The 
results are collected in Tables [1] and [2] in Appendix [Bl 

In Fig. [7] we show (TrU) and (Tr f/~^) as a function of the Langevin stepsize 
for /i = 1 (left) and 3 (right). Statistical fluctuations in the disordered phase 
are larger, even though the Langevin trajectory is twice as long. For the lowest- 
order algorithm stepsize dependence is clearly visible, as expected. The dotted 
lines indicate a linear flt using the data at the four smallest stepsizes. In the 
case of the higher-order algorithm, there appears to be no stepsize dependence 
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Figure 7: Stepsize dependence of (TiUx) (top panes) and (Trt/~^) (bottom 
panes) at = 1 (left) and 3 (right) on a 10^ lattice for /3 = 0.125 and h = 0.02, 
using both the standard lowest-order and the improved algorithm. 
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Figure 8: Stepsize dependence of the real part of (LTrf/j.) and (LTrf/^^) at 
/i = 1 (left) and 3 (right), using both the standard and the improved algorithm. 
Other parameters as in Fig. [71 

visible; the dashed lines indicate the average of the five data points in each case0 
Importantly, we note that the results from both algorithms are consistent in the 
limit e — >■ 0, see also Appendix |Bl 

In order to justify these results, we have computed (LO) (which should be 
equal to zero), where O = TtU,Tt and n. Since these observables are holo- 
morphic, we drop the tilde on the L from now on. Note that (Ln) is not indepen- 
dent, since n is a linear combination of Tr U and Tr U~^. The imaginary parts of 
(LO) are consistent with zero. The stepsize dependence of the real parts is shown 
in Fig. [S] for /i = 1 (left) and 3 (right). Note the different vertical scale: the step- 

^Theoretically, corrections of ©(e'^/^) arc expected [55] . 
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Figure 9: Histograms for KeTrUx (left) and ImTrf/j. (right) at /3 = 0.125, h = 
0.02 and /i = 0, 1,3 on 8^ and 12'^. The vertical lines denote the boundaries in 
SU(3), i.e. without complexification. Note the vertical logarithmic scale. 

size dependence is stronger in the ordered phaseEl For the lowest-order algorithm 
there are again clear finite-stepsize corrections, which vanish in the limit that 
e — 0. In the case of the higher-order algorithm, finite-stepsize corrections are 
much smaller or even absent. We find that (LO) goes to zero in the limit that 
the stepsize is taken to zero. This observation is a necessary requirement for the 
applicability of complex Langevin dynamics. Interestingly, larger finite-stepsize 
corrections correspond to larger deviations of (LO) from zero. It turns out that 
this is also seen when using real Langevin dynamics, e.g. in the phase-quenched 
theory. We conclude therefore that computations of {LO) yield a sensitive test 
to quantify finite-stepsize errors. 

As a final assessment, we discuss the extent to which the complexified field 
space is explored. A sufficiently localized distribution in the imaginary field 
direction is required for the formal justification to hold [3Tl[3l]. In Fig. [9] we 
show histograms for KeTiU (left) and ImTrt/ (right), obtained by binning the 
data sampled during the Langevin process, for three different /i values and two 
lattice volumes at /9 = 0.125. For real dynamics, i.e. when the angles 01^2 are real 
and U G SU(3), Tr U is complex- valued, taking values in a triangular shape with 
corners at 3e^'^'^*/^ (q = 0, 1,2). The corresponding boundaries are shown in the 
figures as vertical dashed lines and the histograms at /i = are contained within 
these boundaries. After complexification, 0i^2 are no longer proper angles and the 
dynamics takes place in the larger group SL(3,C). At nonzero n, we observe that 
the SU(3) boundaries are indeed crossed, as required, but that the distribution 
appears to remain localized (note the vertical logarithmic scale). The tails of 
the distributions are noisy, as they are visited during the Langevin process very 
rarely. There is very little volume dependence. We also note that the histograms 

^Larger stepsize corrections in the ordered phase at larger values of fj, are also seen in the 
Bose gas, where the stepsize is effectively enhanced as e^e [29] , 
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Figure 10: Histogram P(0^), where 0^ = {01,02}, at /3 = 0.125, h = 0.02 and 
fi = 0,1,3 on 8^ and 12^. When = 0, 0^ = 0. The dashed straight hnes are 
P{(j)^) ~ e"''''^ ' with h = 35,45. Note the vertical logarithmic scale. 



at fi = 1 (in the disordered phase) resemble the histograms at /i = 0, while at 
/i = 3 (in the ordered phase) they are significantly different, which is reflected in 
the larger expectation value of (Tr U) . The histograms for Im Tr U are symmetric 
within numerical uncertainty, since (Tr U) is real. 

Finally, in Fig. [TO]we show the histogram for 0^ = {0^, For real Langevin 
dynamics at = 0, 0^ = 0. At flnite /i, nonzero values are generated by the 
complex drift term. We observe that the distribution drops exponentially, over 
many decades, before the signal becomes noisy. The straight dashed lines indicate 
P(0^) ~ e"^'*^ ' with b = 35,45, and are meant to guide the eye. Note that 
the exponential drop is considerably faster than in the U(l) model studied in 
Refs. [3ni3l] , where 6 ~ 2 and complex Langevin dynamics failed. This fast drop 
and localization of the distribution is another requirement for the applicability 
of complex Langevin dynamics|§ 



6 Summary and outlook 

In this paper we revisited the SU(3) spin model, an effective dimensionally re- 
duced Polyakov loop model for QCD in the strong-coupling and heavy-quark 
limit, at nonzero chemical potential. To handle the sign problem we employed 

^An open question is what happens to expectation values of the form (Tr U'') with k large. 
These observables eontain terms of the form e^'^'^^ cos(A:0^) and the presence of the rapidly 
oscillating cosine should be taken into consideration, see also Ref. |34j . 
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complex Langevin dynamics, paying special attention to the justification of the 
method. Using analyticity at small /i^ (Taylor series expansion and smoothness 
in /i^), formal criteria for correctness, and localization of distributions in the 
complexified space, we arrive at the conclusion that complex Langevin dynamics 
is reliable in both the ordered and the disordered phase, including the critical 
region. This should be contrasted with the case of the XY model, where correct 
results were obtained in the ordered phase but neither in the disordered phase 
nor in the transition region [32] . In the XY model this failure was detected by an 
apparent lack of analyticity at small and the presence of very broad, slowly 
decaying distributions (as well as by comparing to results obtained in the world 
line formalism [H]). We can therefore conclude that the assessments employed 
here can be used constructively to rule out or support the applicability of com- 
plex Langevin dynamics|l| We emphasize that these tests are generally applicable 
and not specific to the theory considered here. Besides supporting the results 
of complex Langevin dynamics, we found that the criteria for correctness are 
also relevant for real Langevin dynamics, as they show clear sensitivity to finite- 
stepsize effects. In order to eliminate the leading-order stepsize dependence, we 
have successfully implemented a simple higher-order algorithm and found it to 
remove essentially all stepsize dependence in the observables. 

How can the different behaviour of the XY model and the SU(3) model under 
complex Langevin evolution be understood? One of the features distinguishing 
the two is the presence of a non-trivial Haar measure in the SU(3) case. Pre- 
liminary results indicate that it is this measure which leads to more controlled 
complex Langevin dynamics. We hope to come back to this in the near future |48j . 
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A Higher-order algorithm 

In this appendix, we discuss the higher-order algorithm fl3.12p in some more 
detail. In Ref. [3S] the algorithm is constructed by considering a single update 
step. It was also shown that the set f l3.12p is part of a more general update rule. 
In order to complement the analysis of Ref. [35], we demonstrate the algorithm 
here using a simple linear kernel. We emphasize that the analysis in Ref. [35] is 

''Nevertheless, it will still be interesting to compare with results obtained in the flux formu- 
lation [mill]. 
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for general nonlinear drift term. For notational simplicity, we use a single degree 
of freedom 0. 

The goal is to solve to higher accuracy 



(A.l) 



Consider the linear kernel K = —ucf). The solution of the stochastic process, with 
vanishing initial conditions, is then given by 



n-l 



(1 - eujf ^ 'r/i. 



i=0 



leading to 



lim 



■"n'-t'ni c\ 
n^oo 00 2 — eUJ 00 



(A.2) 



(A.3) 



indicating the linear stepsize dependence (we always assume eu < 1] 
Ref. [35] proposes the following update 

= 0n + ^(^K{(j)n) + k^/ean, 



— Vn 



aK{il)n) + bK{tpn) +Vean 



(A.4) 



where the coefficients a, b, k and / are to be determined, and the noise satisfies 

(A.5) 



~ -1 

I 



and 



Note that 

Straightforward substitution in the case of the linear kernel gives 

0n+l = 0n - e'^^n + -/If]^, 

with 



cu = a;(a + 6) ( 1 — ^ew ) , ?7„ = a„ — (a/c + bl)eooan- 



Noting that 



{Vnr]n') = 2(1 - {ak + bl)eoj + - {ak + blf{eojf ) 5, 



(A.6) 
(A.7) 

(A.8) 
(A.9) 

(A.IO) 
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we find, see Eq. (]A.3p . 

I {flifli) If 1 l + a + h-2{ak + hl) 
n->-oo oj 2 — eu! oj \a + 2[a + o) 



(A.ll) 



We can now determine the coefficients and take 

a + b=l, ak + bl = l. (A.12) 

The resulting expectation value is 

lim = i - + . . .y (A. 13) 

with a remaining correction of O(e^). The two conditions above do not fully 
determine the coefficients a, b, k, I. Ref. [55] argues that a further condition 



ak^ + bf = I (A.14) 

follows from minimizing the local truncation error, which is beyond the scope of 
the linear example discussed here. In the main part of the paper, the coefficients 
are taken as 

a = l, 6 = 1 k = 0, / = ! (A.15) 
which satisfies the constraints above. 



B Tables 

In this Appendix we list the results for the stepsize dependence obtained at 
13 = 0.125 and h = 0.02 on a 10^ lattice, for = 1 (Table [1]) and = 3 (Table E]). 
The total Langevin time is 10,000 for fi = 1 and 5,000 for fj^ = 3, after discarding 
the thermalization stage. 

Every table shows the real part of the three observables (Tr f/^^), (Tr f/^^) and 
(n) in the upper part, and the real part of the criteria (LTiUx), {LTtU~^) and 
(Ln) in the lower part, for both the "lowest-order" and the "improved" algorithm. 
Here lowest-order algorithm refers to the standard discretization fl3.10p . which has 
corrections that are linear in the stepsize; improved algorithm refers to the higher- 
order algorithm f l3.12p of Ref. [35]. In the case of the lowest-order algorithm we 
performed a linear extrapolation, using the values at the four smallest stepsizes. 
Since there is very little stepsize dependence left in the case of the improved 
algorithm, the average is shown. 
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n DDI 






U.UUc/OO J- J- J 


lU Wt-oL 


U.UUU 1 o 




U. 04:00 I J-O J 


U.UUc/OOl 11 I 


r\TTA ay 
Ul Llci 


n nnn^ 




U.OOOO^lO j 


n m n9fii'i i 

U.U lUZU J- J- J 


Q 1 (Tr^Ti "I" n T V I 
diHjUi 1 Lliili 




n 9zL9Ql'97'l 


n '?(sn(Si'9i 

U. OUUU 


U.UIUUUIIO 1 




n nnni 




3584l'90'l 


0103Qri9'l 




00005 

\J • \J\J \J\J ^ 


2413f21) 


3600fl6) 


01047fl0) 




extrapolation 


0.2419(19) 


0.3605(15) 


0.010497(92) 




0.001 


0.2410(27) 


0.3593(21) 


0.01046(13) 


improved 


0.00075 


0.2391(20) 


0.3579(15) 


0.01037(10) 


algorithm 


0.0005 


0.2426(25) 


0.3606(19) 


0.01053(12) 




0.00025 


0.2415(34) 


0.3604(26) 


0.01048(17) 




0.0001 


0.2407(23) 


0.3593(18) 


0.01044(11) 




average 


0.2407(11) 


0.35921(84) 


0.010443(53) 




C 


/ r,Tr TI\ 


/ r,Tr Jl^^\ 






n nni 


U.UUOOl 1^ 1 -L J 




U.UUUlOc/^OO j 




n nnn?^ 


\j .\j\jOO \ j 


U.UUtJZiOlUU I 


n nnm ^oc^o^ 




u.uuuo 






n nnnnfi^ii"?'?'! 


algorithm 


0.00025 


-0.00029(80) 


0.00027(78) 


-0.000018(38) 




0.0001 


0.00094(80) 


0.00109(79) 


0.000043(38) 




0.00005 


0.00056(58) 


0.00068(57) 


0.000026(27) 




0.001 


0.00008(74) 


0.00001(72) 


-0.000014(35) 


improved 


0.00075 


0.00079(60) 


0.00049(59) 


0.000022(28) 


algorithm 


0.0005 


-0.00017(74) 


-0.00036(71) 


-0.000024(35) 




0.00025 


0.0004(12) 


0.0001(12) 


0.000001(58) 




0.0001 


0.00053(69) 


0.00026(68) 


0.000010(33) 



Table 1: Stepsize dependence for /i = 1 (disordered phase), /3 = 0.125 and 
h = 0.02 on a lattice of size 10^. 



References 

[1] 0. Philipsen. larXiv:1111.5370l [hep-ph]. 

[2] P. de Forcrand, PoS LAT2009 (2009) 010 |arXiv:1005.0539l [hep-lat]]. 

[3] S. Gupta, PoSLATTICE 2010 (2010) 007 ja rXiv: 1101. 01 091 [hep-lat]]. 

[4] S. Chandrasekharan and U. -J. Wiese, Phys. Rev. Lett. 83 (1999) 3116 
I arXivicond- mat/9902128^ [cond-mat.stat-mech]]. 

[5] J. Bloch, Phys. Rev. Lett. 107 (2011) 132002 |arXiv:1103.3"467l [hep-lat]]. 



20 





£: 
C 


iTr TJ\ 


\ 1 






n nni 




1 . 1 ouoo \cyo J 


n fi7Q7firi (\\ 


lU WcoL 


n nnn?^ 




1 . 1 OOOO I 1 I 


U.UOUUU114: I 


Ui Lit;! 


n nnn^ 

u.uuuo 


1 7nififi('4n'i 

1. 1 UlUO^^U j 






cliH,Ui 1 Liilii 






1 . 1 4:4: tJU 1 4:1 I 


U.UOOUOllO I 




n nnni 


-L. 1 UtitiZjI OU I 


1 . 1 '4:'4:Zj i ) 


68995('1 5) 




00005 


1 70586 f 42) 


1 74561 f37) 


68352fl7) 




extrapolation 


1.70615(27) 


1.74590(24) 


0.68364(11) 




0.001 


1.70605(40) 


1.74571(36) 


0.68360(16) 


iiiiDrovpd 


00075 


1 70514f27) 


1 74486 f 24) 


68324flll 


algorithm 


0.0005 


1.70629(43) 


1.74598(38) 


0.68370(17) 




0.00025 


1.70597(29) 


1.74571(26) 


0.68357(12) 




0.0001 


1.70596(39) 


1.74576(35) 


0.68356(16) 




average 


1.70576(15) 


1.74549(13) 


0.683485(60) 




£: 
C 


/ r, Tr TJ\ 

\±J ±1 \J j 


/ r,Tr 

\ J-/ li Ly j 






n nni 

U.UUl 


U.U^OUI^ lU j 


U.UOZtJOl^c/ 1 J 






U.UUU 1 tJ 


U.UOlc/OlOU I 


U .U04:00 I OtJ I 


n m 981 f'?4'> 




u.uuuo 


u.uic/^^l^yo J 


u.uzi ly i^yu j 




algorithm 


0.00025 


0.0084(14) 


0.0092(13) 


0.00335(55) 




0.0001 


0.0055(11) 


0.0058(10) 


0.00219(42) 




0.00005 


0.0020(12) 


0.0021(12) 


0.00081(50) 




0.001 


0.0070(11) 


0.0056(11) 


0.00189(43) 


improved 


0.00075 


0.00624(88) 


0.00443(85) 


0.00169(35) 


algorithm 


0.0005 


0.0044(12) 


0.0026(12) 


0.00096(48) 




0.00025 


0.00285(67) 


0.00100(65) 


0.00040(27) 




0.0001 


0.0020(10) 


0.0001(10) 


0.00006(42) 



Table 2: As in Table [H for /i = 3 (ordered phase). 



J. Ambjorn, K. N. Anagnostopoulos, J. Nishimura and J. J. M. Verbaarschot, 
JHEP 0210 (2002) 062 |hep-lat/0208 025|. 

Z. Fodor, S. D. Katz and C. Schmidt, JHEP 0703 (2007) 121 
| hep-lat/0701022] . 

S. Ejiri, Phys. Rev. D 77 (2008) 014508 [arXiv: 0706. 35491 [hep-lat]]. 

K. N. Anagnostopoulos, T. Azuma and J. Nishimura, Phys. Rev. D 83 (2011) 
054504 [arXiv: 1009.45041 [cond-mat.stat-mech]] . 

S. Chandrasekharan, PoS LATTICE2008 (2008) 003 [0810.2419 [hep-lat]]. 



21 



[11] D. Banerjee and S. Chandrasekharan, Phys. Rev. D 81 (2010) 125007 
larXiv: 1001 .36481 [hep-lat]]. 

[12] P. de Forcrand and M. Fromm, Phys. Rev. Lett. 104 (2010) 112005 
[a rXiv:0907.1915l [hep-lat]]. 

[13] J. Langelage, S. Lottini and O. Phihpsen, JHEP 1102 (2011) 057 [Erratum- 
ibid. 1107 (2011) 014] [ arXiv: 10 10.09511 [hep-lat]]. 

[14] Y. D. Mercado, H. G. Evertz and C. Gattringer, Phys. Rev. Lett. 106 (2011) 
222001 [arXiv: 1102.3096 [hep-lat]]. 

[15] C. Gattringer, Nucl. Phys. B 850 (2011) 242 |arXiv: 1104. 25031 [hep-lat]]. 

[16] Y. Delgado, H. G. Evertz, C. Gattringer and D. Goschl, larXiv:1111.0"9T6l 
[hep-lat] . 

[17] W. linger and P. de Forcrand. l arXiv: 111 1.14341 [hep-lat]. 

[18] M. Fromm, J. Langelage, S. Lottini and O. Philipsen. larXivillll. 49531 [hep- 
lat]. 

[19] G. Parisi, Phys. Lett. B 131 (1983) 393. 

[20] J. R. Klauder, Stochastic quantization, in: H. Mitter, C.B. Lang (Eds.), 
Recent Developments in High- Energy Physics, Springer- Verlag, Wien, 1983, 
p. 351; J. Phys. A: Math. Gen. 16 (1983) L317; Phys. Rev. A 29 (1984) 
2036. 

[21] P. H. Damgaard and H. Hiiffel, Phys. Rept. 152 (1987) 227. 

[22] F. Karsch and H. W. Wyld, Phys. Rev. Lett. 55 (1985) 2242. 

[23] N. Bilic, H. Gausterer and S. Sanielevici, Phys. Rev. D 37 (1988) 3684. 

[24] J. Berges and I.-O. Stamatescu, Phys. Rev. Lett. 95 (2005) 202003 
|hep-lat/050 8030|. 

[25] J. Berges, S. Borsanyi, D. Sexty and I. O. Stamatescu, Phys. Rev. D 75 
(2007) 045007 |hep-lat/ 0609058|. 

[26] J. Berges and D. Sexty, Nucl. Phys. B 799 (2008) 306 [0708.0779 [hep-lat]]. 

[27] G. Aarts and I.-O. Stamatescu, JHEP 0809 (2008) 018 [0807.1597 [hep-lat]]. 

[28] G. Aarts, Phys. Rev. Lett. 102 (2009) 131601 [0810.2089 [hep-lat]]. 

[29] G. Aarts, JHEP 0905 (2009) 052 [0902.4686 [hep-lat]]. 



22 



G. Aarts, F. A. James, E. Seller and I. -O. Stamatescu, Phys. Lett. B 687 

(2010) 154 |arXiv:0912.06T7l [hep-lat]]. 

G. Aarts, E. Seller and 1. -O. Stamatescu, Phys. Rev. D 81 (2010) 054508 
[irXlv:0912.3360l [hep-lat]]. 

G. Aarts and F. A. James, JHEP 1008 (2010) 020 [arXlv: 1005.34681 [hep- 
lat]]. 

G. Aarts and K. Splittorff, JHEP 1008 (2010) 017 [arXiv: 1006.03321 [hep- 
lat]]. 

G. Aarts, F. A. James, E. Seller and I. -O. Stamatescu, Eur. Phys. J. C 71 

(2011) 1756 [arXlv:1101.3270l [hep-lat]]. 

Chien-Cheng Chang, Math. Comp. 49 180 (1987) 523-542. 

G. Aarts, F. A. James, E. Seller and I. -O. Stamatescu, PoSLATTICE 2011 
(2011) 197 [arXiviim5749. [hep-lat]]. 

T. D. Cohen, Phys. Rev. Lett. 91, 222001 (2003) |hep-ph/0307089] . 

J. C. Osborn, K. Splittorff and J. J. M. Verbaarschot, Phys. Rev. Lett. 94, 
202001 (2005) |hep-th/0 5012l0l. 

P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential 
Equations, Springer Verlag (1992). 

I. T. Drummond, S. Duane and R. R. Horgan, Nucl. Phys. B 220 (1983) 
119. 

S. M. Catterall, L T. Drummond and R. R. Horgan, Phys. Lett. B 254 
(1991) 177. 

J. R. Klauder and W. P. Petersen, J. Stat. Phys. 39 (1985) 53. 

J. Ambjorn and S. K. Yang, Phys. Lett. B 165 (1985) 140. 

J. Ambjorn, M. Flensburg and C. Peterson, Nucl. Phys. B 275 (1986) 375. 

K. Okano, L. Schulke and B. Zheng, Prog. Theor. Phys. Suppl. Ill (1993) 
313. 



M. -P. Lombardo, Nucl. Phys. Proc. Suppl. 83 (2000) 375 |hep-lat/9908006] . 



P. de Forcrand and O. Philipsen, Nucl. Phys. B 642 (2002) 290 
| hep-lat/0205016] . 

[48] G. Aarts, F. A. James, E. Seller, D. Sexty and I. -O. Stamatescu, in prepa- 
ration. 



23 



